Quantum molecular dynamic simulations of warm dense carbon monoxide 
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Using quantum molecular dynamic simulations, we have studied the thermophysical properties of 
warm dense carbon monoxide under extreme conditions. The principal Hugoniot, which is derived 
from the equation of state, shows excellent agreement with available experimental data up to 67 
GPa. The chemical decomposition of carbon monoxide has been predicted at 8 GPa by means of 
pair correlation function. Based on Kubo-Greenwood formula, the dc electrical conductivity and the 
optical reflectivity are determined, and the nonmetal-metal transition for shock compressed carbon 
monoxide is observed around 43 GPa. 
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I. INTRODUCTION 

With rapid development of planetary science and ex- 
plosion technology, pressure induced responses of ma- 
terials under extreme conditions, where the combina- 
tion of high temperatures and high pressures defines 
the warm dense regime, have recently attracted exten- 
sive attentioE^ii"— . Experimental techniques, e.g. gas 
guni, chemical explosives^, pulsed power— li, have been 
adopted in dynamic compressions, where pressures could 
reach megabar region. Due to the enormous progress 
in computational capacity, quantum molecular dynam- 
ics (QMD)*'^, provide powerful tools to study warm 
dense matter theoretically, where dissociations, recom- 
binations, and ionizations characterize the high pressure 
behaviors. 

As one of the major components in the reacted high 
explosives, the pressure introduced transitions of carbon 
monoxide (CO) are of crucial interest in exploring deto- 
nation products^"—. The equation of states (EOS) and 
the relative quantities, for instance the Hugoniot points, 
sound velocity and particle velocity, are important pa- 
rameters for materials under extreme conditions. Fur- 
thermore, electrical and optical properties such as, e.g., 
dc conductivity (adc) and optical reflectance are closely 
related to the dielectric function, which can be evalu- 
ated through dynamic conductivity. The above men- 
tioned quantities are of crucial interest in characterizing 
the unique behavior of warm dense CO. The Hugoniot 
points up to 70 GPa have been detected by a two-stage 
light-gas gun^, and three theoretical models were used to 
analyze the EOS and chemical decomposition. CO was 
reported to decomposes into condensed carbon and fluid 
CO2 at high pressures by pair potential calculations^^, 
and the dissociation rates were then determined-i^. The 
principal Hugoniot and decomposition of CO have been 
investigated by means of QMD method^. However, wide 
range EOS, shock products, electrical and optical prop- 
erties of CO under extreme conditions are still yet to be 
presented and understood. 

In the present work, we have performed QMD simu- 
lations to study the thermophysical properties of warm 



dense CO. The decomposition and recombination of CO 
along the principal Hugoniot are analyzed by pair corre- 
lation functions (PCF). The dynamic conductivity is cal- 
culated by Kubo-Greenwood formula, from which the dc 
conductivity, the dielectric function and optical reflectiv- 
ity can be extracted. The rest of the paper is arranged 
as follows. In section II, the computational method is 
briefly described. We present and discuss our calculated 
results, some of which are made comparison with avail- 
able experimental and theoretical results in section HI. 
Finally, we close our paper with a summary of our main 
results. 



II. COMPUTATIONAL METHOD 

In this study, the molecular dynamics trajectories are 
calculated by employing the Vienna ab initio simulation 
package (VASP) plane-wave pseudopotential code devel- 
oped at the Technical University of Viennai^!^. The 
calculations are performed by a series of volume-fixed su- 
percell with an invariable number of atoms. Employing 
the Born-Oppenheimer approximation, the electrons are 
fully quantum-mechanical treated by solving the Kohn- 
Sham equations for the orbitals and energies within 
a plane-wave, finite-temperature DFT formulation^-, 
where the electronic states are populated according to 
the Fermi-Dirac distribution at electron temperature 
Te. The Perdew-Wang 91 (PW91) generalized gradi- 
ent approximation— and the projector augmented wave 
(PAW) potential of Blochli^ are employed to describe 
the exchange-correlation energy and the electron-ion in- 
teraction, respectively. The isokinetic ensemble (NVT) 
is employed for the ions, where the ionic temperature Ti 
is controlled by Nose thermostat^- and the system is con- 
trolled in local thermodynamical equilibrium by setting 
the electron temperature Te and the ion temperature Ti 
to be equal. 

The supercell in our calculation contains 128 atoms (64 
CO molecules) with periodic boundary condition. The 
plane-wave cutoff energy is selected to be 550 eV, which 
is tested to be good convergence for both total energy 
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TABLE I: Hugoniot pressure (P) and temperature (T) points 
derived from QMD simulations at a series of densities (p). 
The corresponding particle velocity (up), and shock velocity 
(its) are also displayed. 



p 


P 


T 


Up 


Us 


(g/cm^) 


(GPa) 


(K) 


(km/s) 


(km/s) 


1.45 


5.32 


603 


1.70 


3.83 


1.78 


8.57 


1614 


2.40 


4.39 


2.25 


20.20 


4545 


4.00 


6.23 


2.52 


31.77 


6600 


5.17 


7.60 


2.70 


43.27 


8697 


6.12 


8.74 


2.90 


53.30 


9175 


6.90 


9.56 


3.10 


67.27 


11000 


7.84 


10.61 



and pressure. In molecular dynamic simulation, only F 
point of the Brillouin zone is included, while, 4x4x4 
Monkhorst-Pack scheme^S grid points are used in the 
electronic structure calculations. The selected densities 
range from 0.807 to 3.10 g/cm^ with temperatures from 
77 to 11000 K, which highlight the principal shock Hugo- 
niot region. All the dynamic simulations last 3-5 ps with 
time steps of 0.5-1.0 fs according to different conditions. 
The EOS data and pair correlation functions are obtained 
by averaging over all the particles and the final 1-2 ps 
simulations. 



III. RESULTS AND DISCUSSION 
A. Equation of state and pair correlation function 

A precise description of EOS plays an important role in 
accurately calculating the electrical and optical proper- 
ties. The EOS is examined theoretically along the Hugo- 
niot, which can be derived from conservation of mass, 
momentum, and energy for a isolated system compressed 
by a pusher at a constant velocity. Rankine-Hugoniot 
equation describes the shock adiabat through a relation 
between the initial and final internal energies, pressures 
and volumes as follows^: 

(E, ^Eo) + ^ {Vi - Vo) (Po + Pi) = 0, (1) 

where V and P are the internal energy, volume, and 
pressure, and the subscripts and 1 denote the initial and 
shocked states, respectively. In our canonical ensemble 
calculations, the internal energy E equals to the sum of 
the total energy from the finite-temperature DFT calcu- 
lation and zero-point energy. The pressure P consists of 
the electronic Pe and ionic Pi components. The pressure 
can be obtained by P = Pg + puksT, where /?„ is the ion 
number density, and fcs is the Boltzman constant. 

As the starting point along the Hugoniot, the ini- 
tial density is po=0.807 g/cm'^ at a temperature of 77 
K^, with the relative internal energy Eq — —14.84 
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FIG. 1: (Color online) Principal Hugoniot (the main panel) 
and the {us, Up) diagram (the inset) of CO. For comparison, 
the present QMD results, together with the previous exper- 
imental data (Ref. 4) and other theoretical results (Ref. 4 
and Ref. 12) are all provided. 



eV/molecule. The initial pressure can be neglected com- 
pared to the high pressure of the final state. The 
Hugoniot points are determined in the following way. 
For a given density p, a series of simulations are ex- 
ecuted for different temperature T. {Ei — Eq) and 
(Vb — Vi)(Po + Pi)/2 as a function of temperature are 
then fitted to polynomial expansions of T. The principal 
Hugoniot temperature Ti and pressure Pi can be deter- 
mined by solving Eq. (1). The particle velocity Up and 
shock velocity Us are then determined from the other two 
Rankine-Hugoniot equationsSi^ Vi — Vb [1 — {up/us)] 
and Pi — Pq = poUgUp. The principal Hugoniot points of 
CO derived from Eq. (1) as well as the particle velocity 
Up and shock velocity Ug are listed in Table I. 

The simulated principal Hugoniot points for CO are 
shown in Fig. 1, where previous experimental and theo- 
retical data are also provided for comparison. As seen in 
Fig. 1, our simulated results accord excellently with the 
experimental results* and the theoretical results^. The 
(up, Us) diagram also shows good agreement with exper- 
imental results (the inset of Fig. 1). 

To examine the structural transitions of CO under 
shock compressions, the PCFs [g{r)= ^ tv(R,R+dR) ] ' 
which represent the possibility of finding a particle at 
a distance r from a reference atom, are calculated. The 
PCFs and the corresponding atomic structures along the 
principal Hugoniot of CO, are presented in Fig. 2 and 
Fig. 3, respectively. At the pressure of 5 GPa (p=1.45 
g/cm3 and T=603 K), the peak of C-0 PCF gc-o(7-) 
locates around 1.13 A, which corresponds to the equi- 
librium internuclear distance of the C-0 bond in CO 
molecule^^ [Fig. 2(a)], indicating that CO molecules re- 
main their ideal molecular fluid state without dissoci- 




FIG. 2: (Color online) Pair-correlation functions for C-C 
(black line), O-O (red line), C-O (blue line) along the prin- 
cipal Hugoniot of CO. 



ation. From the corresponding atomic structure [Fig. 
3(a)], it can be clearly seen that the supercell is filled 
with CO molecules. With the increase of pressure along 
the Hugoniot, the main peak of PCF gc-o(^) becomes 
reduced in amplitude and gets broadened due to the ther- 
mal excitation, as shown in Fig. 2(b). Meanwhile, the 
gc-c('') appears a peak at r=1.47A, which lies between 
typical sp^ and sp^ hybrid C-C bond length. From the 
atomic structure of CO at p=2.2b g/cm^ (P=20.2 Gpa, 
T=4545 K), CO2 molecules and the small clusters with 
carbon backbones can be clearly observed [Fig. 3(b)]. 
The structural transitions of CO may bring with the in- 
crease of the adc , which will be discussed in the following 
section. As the difference between the length of C-O 
bond in CO molecule (1.13A) and that of CO2 molecule 
(1.16A)23 is very small, the appearances of their PCFs 
are similar. Despite the fact that the formation of CO2 
increases the peak value of gc-o(''), CO is further dis- 
sociated leading to a decrease in the peak value. Under 
the competition of the above two effects, the amplitude 
of PCF gc-o('^) still keeps reduced as the pressure in- 
creases. At higher pressures, CO2 molecules get dissoci- 
ated, and a mixture of CO, atomic carbon, and some car- 
bon backbones clusters forms, which results in a further 
decrease of the peak value of PCF gc-o('")- In particu- 
lar, diatomic oxygen [the peak of PCF go-o(^) is around 
1.23 A as shown in Fig. 2(d)] and monoatomic oxygen are 
formed in the shocked system [Fig. 3(d)], which lead to 
the nonmetal-metal transition. Overall, the systematic 
behavior of the dissociation and recombination of fluid 
CO under shock pressure is described. At low pressures 
(P<20 GPa), CO2 as well as some carbon backbone clus- 
ters forms. With increasing the pressure (P>43GPa), 
CO2 dissociates and a mixture of CO, atomic carbon di- 




FIG. 3: (Color online) Atomic structures at different Hugo- 
niot points. The carbon and oxygen atoms are denoted by 
gray and red balls, respectively, and the relative isosurface of 
charge density (blue regimes) are plotted. (a)p=1.45 g/cm'^, 
r=603 K; (b)p=2.25 g/cm^ r=4545 K; (c)p=2.52 g/cm^ 
r=6600 K; (d)p=2.70 g/cm^ r=8697 K. 



atomic oxygen monoatomic oxygen and carbon backbone 
clusters forms. 



B. Dynamic conductivity under shock 
compressions 

The dissociation of molecules under dynamic compres- 
sion are closely related to the nonmetal-metal transition 
properties. The dynamic conductivity is derived from the 
Kubo-Greenwood formula as foUow o^'^i^^ : 



fJi (w) 



2tt 



N N 3 

^^k)5]5]^[/(e„k)-/(e„k)] 

Ic J — 1 i— la— 1 



(*j,k|V„|1'^,k)|^(5(ej,k - ei,k - hoj), (2) 



where ^'i^k is Kohn-Sham eigenstate, with corresponding 
eigenvalue ei_k, and occupation number /(ei,k). fl is the 
volume of the supercell. u (k) is the K-point weighting 
factor. The i and j summations range over N discrete 
bands included in the calculation, and a over the three 
spatial directions. 

In order to calculate the conductivity, thirty atomic 
configurations are selected at equilibrium along the prin- 
cipal Hugoniot. The behavior of the conductivity cti ver- 
sus energy at different Hugoniot points are calculated 
(not shown). With the increase of density and temper- 
ature along the principal Hugoniot, the main peak in- 
creases in amplitude and gets broadened. Such change 
of the peak indicates the increase of another impor- 
tant physical quantity, dc conductivity, which follows the 
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FIG. 4: The dc conductivity dependence on pressure is ex- 
tracted from the dynamic electrical conductivity. 
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FIG. 5: (Color online) Optical reflectivity of shocked CO 
for wavelengths of 404, 808 and f064 nm along the princi- 
pal Hugoniot. 



static limit tJdc — limai (uj). The dc conductivity is ex- 

tracted as shown in Fig. 4. For the pressures below 8 
GPa, (Tdc can be treated as small as negligible, which 
shows the insulating nature of CO at low pressures. At 
higher pressures, the dc conductivity adc of CO jumps 
three orders of magnitude from 8 to 67 GPa (603 to 
11000 K). The nonmetal-metal transition occurs around 
43 Gpa, which is accompanied by the gradual dissocia- 
tion of the molecular fluid. Such nonmetal-metal transi- 
tion under shock pressure has also been reported in fluid 
COz^ and Oz^. 

The imaginary part of conductivity a2 can be obtained 
via the Kramers-Kronig relation as 



0-2 (w) 



-P 



<Jl(v)u! 
(zy2 -^2) 



dv, 



(3) 



where P is the principal value of the integral and v is 
frequency. The real part and the imaginary part 

e2(w) of dielectric functions can be derived from the two 
parts of the conductivity. 



and 



ei(w) = 1 ^-cr2(w) 



£2(1^) = cri(a;) 

€qUJ 



(4) 



(5) 



The square of the index of refraction, which contains 
the real part n(a;) and the imaginary part k{ui), is equal 
to the dielectric function. Then the index of refraction 
can be obtained by 



and 



k{u)^\l-Mu:)\-e,{uj)]. 



(6) 



(7) 



The index of refraction is useful for evaluating optical 
properties, such as the reflectivity and absorption coeffi- 
cient. Emissivity of spectrum is an important quantity, 
which can be used to determine temperature in exper- 
iments. The optical reflectivity and emissivity of spec- 
trum are closely related. We calculated the reflectivity 
of CO under shock pressure by the following relation 



r{uj) = 



[I - n{uj)Y + k{ujf 
[I + n{Lo)Y + kiuf 



(8) 



The optical reflectivities of shocked CO for wavelengths 
of 404, 808 and 1064 nm along the principal Hugoniot are 
shown in Fig. 5. The change induced by the pressure can 
be clearly seen. The optical reflectance increases from 
0.12 to 0.46-0.60 with the pressure increasing from 8 to 
67 GPa, which can be interpreted as a gradual transition 
from nonmetal state to metal state. Similar behavior 
in the optical reflectance of liquid deuterium has been 
previously observed in experimental. Our predictions of 
CO need be tested in the future experiments. 



IV. CONCLUSIONS 

In summary, we have performed QMD simulations to 
study the thermophysical and optical properties of shock- 
compressed carbon monoxide along the principal Hugo- 
niot. The principal Hugoniot curve has been derived 
from the EOS, which shows excellent agreement with ex- 
perimental and theoretical results. The dissociation of 
CO under shock compressions is described by the pair 
correlation functions. At low pressures (P<20 GPa), 
CO2 and large carbon backbone clusters form. Above 
43 GPa, CO2 molecules dissociate and a mixture of 
CO, atomic carbon and carbon backbones clusters form. 
Based on the calculation of the dynamic conductivity, the 
nonmetal-metal transition and the change in optical re- 
flectivity are observed. In particular, the nonmetal-metal 
transition occurs around 43 GPa. 
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